{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<a target=\"_blank\" href=\"https://colab.research.google.com/github/opensim-org/opensim-core/blob/main/Bindings/Python/tutorials/Tutorial%207%20-%20Introduction%20to%20OpenSim%20Moco.ipynb\">\n",
        "  <img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/>\n",
        "</a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YS-0eq5Y6BjU"
      },
      "source": [
        "# 7. Introduction to OpenSim Moco"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZZClpIjbitD4"
      },
      "source": [
        "\n",
        "## 7.1. Objectives"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AKANMMa6iycq"
      },
      "source": [
        "**Introduction**\n",
        "\n",
        "[OpenSim Moco](https://opensim-org.github.io/opensim-moco-site/) is a software package for creating musculoskeletal simulations using trajectory optimization [[1]](https://doi.org/10.1371/journal.pcbi.1008493). Users define problems by providing an OpenSim model, cost function terms, and additional problem constraints via an easy-to-use scripting interface. Moco rapidly solves trajectory optimization problems via the state-of-the-art direct collocation method.\n",
        "\n",
        "**Purpose**\n",
        "\n",
        "The purpose of this tutorial is to demonstrate how to use Moco to build, solve, and analyze a simple trajectory optimization problem using the [OpenSim conda package](https://anaconda.org/opensim-org/opensim) and [Jupyter notebooks](https://jupyter.org/). In this tutorial you will:\n",
        "\n",
        "*   Set up a new conda environment and install OpenSim with Moco.\n",
        "*   Learn the basics of Moco's `MocoStudy` interface for defining optimal control problems.\n",
        "*   Solve Moco's simplest example, the \"sliding mass\" problem.\n",
        "*   Plot the results using `matplotlib`."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Sliding Mass Problem**\n",
        "\n",
        "You will use Moco to find the optimal trajectory of a point mass starting at rest, translating one meter, and ending at rest, in minimum time.\n",
        "\n",
        "\n",
        "![sliding_mass](https://drive.google.com/uc?export=download&id=17YTG2WqP6Hzg83OwGyOjqXmyqIWyl34s)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Local Installation**\n",
        "\n",
        "You may also run this notebook using a locally installed OpenSim `conda` package. Visit [the OpenSim python scripting Confluence page](https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53085346/Scripting+in+Python) for more information on how to install OpenSim in a local Python environment. You will also need to download the resources from Google Drive to the same directory as this notebook. Then, you may skip to Section 7.3 below."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kPnGlZIT8oiM"
      },
      "source": [
        "## 7.2. Install `conda` and OpenSim"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3kxYrs4Hi6Nl"
      },
      "source": [
        "In the following cell, you will use [condacolab](https://github.com/conda-incubator/condacolab) to install [Miniconda](https://docs.conda.io/en/latest/miniconda.html). If you want to install a different conda distribution (e.g., Anaconda, Mambaforge) check [condacolab's documentation](https://github.com/conda-incubator/condacolab#usage)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "5-sX9L3f6Asq",
        "outputId": "07c714c8-7744-4bf4-fa0f-f1c66d2fcec9"
      },
      "outputs": [],
      "source": [
        "!pip install -q condacolab\n",
        "import condacolab\n",
        "condacolab.install()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LdmpfqC8Dt0x"
      },
      "source": [
        "Don't worry if after executing the previous cell you get an error saying that your session has failed. This happens because `condacolab` needs to restart the session for the changes to take effect. Therefore, you have to execute the previous cell before executing any other cell.\n",
        "\n",
        "Now, execute the following cell to install the OpenSim `conda` package."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "H9WE2h3fD61J",
        "outputId": "2795df96-b021-40eb-ba47-39a105cc2ac6"
      },
      "outputs": [],
      "source": [
        "!conda install opensim-org::opensim"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 7.3. Import the OpenSim module"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BV3ofVUoEkEi"
      },
      "source": [
        "The following cell imports the OpenSim module and, if successful, shows the current version and build date."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "import opensim as osim\n",
        "osim.GetVersionAndDate()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 7.4. Create the sliding mass model\n",
        "\n",
        "First, let's create an OpenSim `Model` representing a \"sliding mass\": a point mass model with a 1-DOF `SliderJoint` actuated by a `CoordinateActuator`. Moco will use this model to enforce the dynamics defect constraints in the direct collocation problem we will solve below."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "uCT1o9zzF-Nn"
      },
      "outputs": [],
      "source": [
        "model = osim.Model()\n",
        "model.setName('sliding_mass')\n",
        "model.set_gravity(osim.Vec3(0, 0, 0))\n",
        "body = osim.Body('body', 2.0, osim.Vec3(0), osim.Inertia(0))\n",
        "model.addComponent(body)\n",
        "\n",
        "# Allows translation along x.\n",
        "joint = osim.SliderJoint('slider', model.getGround(), body)\n",
        "coord = joint.updCoordinate()\n",
        "coord.setName('position')\n",
        "model.addComponent(joint)\n",
        "\n",
        "actu = osim.CoordinateActuator()\n",
        "actu.setCoordinate(coord)\n",
        "actu.setName('actuator')\n",
        "actu.setOptimalForce(1)\n",
        "model.addComponent(actu)\n",
        "\n",
        "body.attachGeometry(osim.Sphere(0.05))\n",
        "\n",
        "model.finalizeConnections()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 7.5. Create a `MocoStudy`\n",
        "\n",
        "`MocoStudy` is the most flexible interface for defining optimal control problems in Moco. It encapsulates both a `MocoProblem`, the optimal control problem definition, and a `MocoSolver`, the module which converts the optimal control problem to a nonlinear program using direct collocation.\n",
        "\n",
        "Run the following code cell to create a new `MocoStudy`, access the underlying `MocoProblem`, and assign the problem the `Model` we created above."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Create MocoStudy.\n",
        "# ================\n",
        "study = osim.MocoStudy()\n",
        "study.setName('sliding_mass')\n",
        "\n",
        "# Define the optimal control problem.\n",
        "# ===================================\n",
        "problem = study.updProblem()\n",
        "\n",
        "# Model (dynamics).\n",
        "# -----------------\n",
        "problem.setModel(model)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 7.6. Define the state bounds and objective function\n",
        "\n",
        "The sliding mass problem is a boundary value problem (BVP), so we need to set some bounds for the state variables, e.g., the mass position and speed. Run the following cell to set the initial position to 0, the final position to 1, and the initial and final speeds to zero. \n",
        "\n",
        "This cell also configures a minimum time objective function using `MocoFinalTimeGoal`. Since we are optimizing time, the final time bounds are set to the interval `[0, 5]` seconds."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Bounds.\n",
        "# -------\n",
        "# Initial time must be 0, final time can be within [0, 5].\n",
        "problem.setTimeBounds(osim.MocoInitialBounds(0.), osim.MocoFinalBounds(0., 5.))\n",
        "\n",
        "# Position must be within [-5, 5] throughout the motion.\n",
        "# Initial position must be 0, final position must be 1.\n",
        "problem.setStateInfo('/slider/position/value', osim.MocoBounds(-5, 5),\n",
        "                     osim.MocoInitialBounds(0), osim.MocoFinalBounds(1))\n",
        "# Speed must be within [-50, 50] throughout the motion.\n",
        "# Initial and final speed must be 0. Use compact syntax.\n",
        "problem.setStateInfo('/slider/position/speed', [-50, 50], [0], [0])\n",
        "\n",
        "# Applied force must be between -50 and 50.\n",
        "problem.setControlInfo('/actuator', osim.MocoBounds(-50, 50))\n",
        "\n",
        "# Cost.\n",
        "# -----\n",
        "problem.addGoal(osim.MocoFinalTimeGoal())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 7.7. Solve the direct collocation problem\n",
        "\n",
        "Next, we will use `MocoCasADiSolver` to transcribe our optimal control problem to a nonlinear program (NLP) using direct collocation. This `MocoSolver` uses CasADi [[2]](https://doi.org/10.1007/s12532-018-0139-4), an open-source software framework for numerical optimization. Internally, Moco uses CasADi to convert your `MocoProblem` to a direct collocation problem written using CasADi's symbolic and numerical types.\n",
        "\n",
        "Run the following cell to intialize the `MocoCasADiSolver`, solve the the optimal control problem, and write the solution to a file. Note that if you want to re-run this cell, you will need to reset the Jupyter notebook to clear workspace variables since `study.initCasADiSolver()` cannot be called after a `MocoSolver` has already been initialized."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Configure the solver.\n",
        "# =====================\n",
        "solver = study.initCasADiSolver()\n",
        "solver.set_num_mesh_intervals(100)\n",
        "\n",
        "# Now that we've finished setting up the study, print it to a file.\n",
        "study.printToXML('sliding_mass.omoco')\n",
        "\n",
        "# Solve the problem.\n",
        "# ==================\n",
        "solution = study.solve()\n",
        "\n",
        "solution.write('sliding_mass_solution.sto')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bMYcGGuaetwS"
      },
      "source": [
        "You should now have `sliding_mass.omoco` and `sliding_mass_solution.sto` in your files."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 7.8. Plot the solution\n",
        "\n",
        "Finally, let's plot the solution. After solving the problem above, `MocoStudy` returned a `MocoSolution` containing the state and control trajectories which we can access directly in Python using `getStateMat()` and `getControlMat()`, respectively."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "\n",
        "time = solution.getTimeMat()\n",
        "position = solution.getStateMat('/slider/position/value')\n",
        "speed = solution.getStateMat('/slider/position/speed')\n",
        "actuator = solution.getControlMat('/actuator')\n",
        "\n",
        "# Create six subplots, with 2 rows and 3 columns.\n",
        "fig, axs = plt.subplots(1, 3, figsize=(12, 5))\n",
        "fig.suptitle('Sliding Mass Solution')\n",
        "\n",
        "# Plot the knee angles on the first subplot.\n",
        "axs[0].plot(time, position, label='position', lw=4)\n",
        "axs[0].set_ylabel('position (m)')\n",
        "axs[0].set_xlabel('time (s)')\n",
        "\n",
        "axs[1].plot(time, speed, label='speed', lw=4)\n",
        "axs[1].set_ylabel('speed (m/s)')\n",
        "axs[1].set_xlabel('time (s)')\n",
        "\n",
        "axs[2].plot(time, actuator, label='actuator', lw=4)\n",
        "axs[2].set_ylabel('actuator force (N)')\n",
        "axs[2].set_xlabel('time (s)')\n",
        "\n",
        "for ax in axs:\n",
        "    ax.set_xlim(0, time[-1])\n",
        "    ax.grid()\n",
        "\n",
        "plt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Questions**\n",
        "\n",
        "*  Describe the state and control trajectories returned after solving the sliding mass problem. Given the optimal control problem defined above, does this solution make sense to you? Explain. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "J5mhS9zNX9ff"
      },
      "source": [
        "## 7.9. Conclusion"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_M04jPJLuiUz"
      },
      "source": [
        "In this tutorial you used [condacolab](https://github.com/conda-incubator/condacolab) to install [Miniconda](https://docs.conda.io/en/latest/miniconda.html) in [Google Colab](https://colab.research.google.com/?hl=en). Then, you installed Conda [OpenSim](https://opensim.stanford.edu/) with Moco into the Miniconda environment. You used the OpenSim API to build a \"sliding mass\" model, and then used Moco to define an optimal control problem using the `MocoStudy` interface. Finally, you plotted the results using `matplotlib`."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lLWgd8Noun6Y"
      },
      "source": [
        "## 7.10. Useful Links"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "m8M3dlT5xQcD"
      },
      "source": [
        "> **OpenSim Website:** https://opensim.stanford.edu/\n",
        ">\n",
        "> **OpenSim Moco Website:** https://opensim-org.github.io/opensim-moco-site/\n",
        ">\n",
        "> **OpenSim Python Scripting:** https://opensimconfluence.atlassian.net/wiki/spaces/OpenSim/pages/53085346/Scripting+in+Python\n",
        ">\n",
        "> **OpenSim API Documentation:** https://simtk.org/api_docs/opensim/api_docs/\n",
        ">\n",
        "> **OpenSim Creator Website:** https://opensimcreator.com/\n",
        ">\n",
        "> **SimTK Website:** https://simtk.org/projects/opensim\n",
        ">\n",
        "> **Biomechanics of Movement Course Videos:** https://www.youtube.com/channel/UCDNGy0KKNLQ-ztcL5h2Z6zA"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rFR2Rf9tk_9A"
      },
      "source": [
        "## 7.11. Acknowledgments\n",
        "\n",
        "Thanks to [OpenSimColab](https://simtk.org/projects/opencolab) project [[3]](https://doi.org/10.1080/10255842.2022.2104607) for creating the first OpenSim Conda package."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dxMjdukAnH2N"
      },
      "source": [
        "## 7.12. References"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WX2A_iaYulV1"
      },
      "source": [
        "> [1] Dembia CL, Bianco NA, Falisse A, Hicks JL, Delp SL (2020) **OpenSim Moco: Musculoskeletal optimal control.** *PLOS Computational Biology* 16(12): e1008493. https://doi.org/10.1371/journal.pcbi.1008493\n",
        ">\n",
        "> [2] Andersson JAE, Gillis J, Horn G, Rawlings JB, Diehl M (2019) **CasADi – A software framework for nonlinear optimization\n",
        "and optimal control.** *Mathematical Programming Computation* 11(1): 1-36. https://doi.org/10.1007/s12532-018-0139-4\n",
        ">\n",
        "> [3] Mokhtarzadeh, H., Jiang, F., Zhao, S., & Malekipour, F. (2022). **OpenColab project: OpenSim in Google colaboratory to explore biomechanics on the web.** *Computer Methods in Biomechanics and Biomedical Engineering*, 1–9. https://doi.org/10.1080/10255842.2022.2104607"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "opensim",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.11.10"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}
